Link to notebook
Link to github repo.
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
[30m── [1mAttaching packages[22m ────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──[39m
[30m[32m✓[30m [34mggplot2[30m 3.3.2 [32m✓[30m [34mpurrr [30m 0.3.4
[32m✓[30m [34mtibble [30m 3.0.4 [32m✓[30m [34mdplyr [30m 1.0.2
[32m✓[30m [34mtidyr [30m 1.1.2 [32m✓[30m [34mstringr[30m 1.4.0
[32m✓[30m [34mreadr [30m 1.4.0 [32m✓[30m [34mforcats[30m 0.5.0[39m
[30m── [1mConflicts[22m ───────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
[31mx[30m [34mdplyr[30m::[32mfilter()[30m masks [34mstats[30m::filter()
[31mx[30m [34mdplyr[30m::[32mlag()[30m masks [34mstats[30m::lag()[39m
library(readxl)
library(phyloseq)
#library(dada2)
library(Biostrings)
Loading required package: BiocGenerics
Loading required package: parallel
Attaching package: ‘BiocGenerics’
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport, clusterMap,
parApply, parCapply, parLapply, parLapplyLB, parRapply, parSapply, parSapplyLB
The following objects are masked from ‘package:dplyr’:
combine, intersect, setdiff, union
The following objects are masked from ‘package:stats’:
IQR, mad, sd, var, xtabs
The following objects are masked from ‘package:base’:
anyDuplicated, append, as.data.frame, basename, cbind, colnames, dirname, do.call,
duplicated, eval, evalq, Filter, Find, get, grep, grepl, intersect, is.unsorted, lapply,
Map, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply, union, unique, unsplit,
which, which.max, which.min
Loading required package: S4Vectors
Loading required package: stats4
Attaching package: ‘S4Vectors’
The following objects are masked from ‘package:dplyr’:
first, rename
The following object is masked from ‘package:tidyr’:
expand
The following object is masked from ‘package:base’:
expand.grid
Loading required package: IRanges
Attaching package: ‘IRanges’
The following object is masked from ‘package:phyloseq’:
distance
The following objects are masked from ‘package:dplyr’:
collapse, desc, slice
The following object is masked from ‘package:purrr’:
reduce
Loading required package: XVector
Attaching package: ‘XVector’
The following object is masked from ‘package:purrr’:
compact
Attaching package: ‘Biostrings’
The following object is masked from ‘package:base’:
strsplit
#library(DECIPHER)
library(phangorn)
Loading required package: ape
Attaching package: ‘ape’
The following object is masked from ‘package:Biostrings’:
complement
library(readr)
library(seqinr)
Attaching package: ‘seqinr’
The following objects are masked from ‘package:ape’:
as.alignment, consensus
The following object is masked from ‘package:Biostrings’:
translate
The following object is masked from ‘package:dplyr’:
count
library(decontam)
library(ape)
library(vegan)
Loading required package: permute
Attaching package: ‘permute’
The following object is masked from ‘package:seqinr’:
getType
Loading required package: lattice
This is vegan 2.5-7
Attaching package: ‘vegan’
The following objects are masked from ‘package:phangorn’:
diversity, treedist
#library(philr)
library(RColorBrewer)
library(microbiome)
microbiome R package (microbiome.github.com)
Copyright (C) 2011-2020 Leo Lahti,
Sudarshan Shetty et al. <microbiome.github.io>
Attaching package: ‘microbiome’
The following object is masked from ‘package:vegan’:
diversity
The following object is masked from ‘package:phangorn’:
diversity
The following object is masked from ‘package:Biostrings’:
coverage
The following objects are masked from ‘package:IRanges’:
coverage, transform
The following object is masked from ‘package:S4Vectors’:
transform
The following object is masked from ‘package:ggplot2’:
alpha
The following object is masked from ‘package:base’:
transform
#library(DESeq2)
library(compositions);
Welcome to compositions, a package for compositional data analysis.
Find an intro with "? compositions"
Attaching package: ‘compositions’
The following object is masked from ‘package:seqinr’:
alr
The following object is masked from ‘package:ape’:
balance
The following objects are masked from ‘package:IRanges’:
cor, cov, var
The following objects are masked from ‘package:S4Vectors’:
cor, cov, var
The following objects are masked from ‘package:BiocGenerics’:
normalize, var
The following objects are masked from ‘package:stats’:
cor, cov, dist, var
The following objects are masked from ‘package:base’:
%*%, norm, scale, scale.default
library(cowplot)
library(plotly)
Attaching package: ‘plotly’
The following object is masked from ‘package:XVector’:
slice
The following object is masked from ‘package:IRanges’:
slice
The following object is masked from ‘package:S4Vectors’:
rename
The following object is masked from ‘package:ggplot2’:
last_plot
The following object is masked from ‘package:stats’:
filter
The following object is masked from ‘package:graphics’:
layout
library(htmlwidgets)
library(withr)
library(lubridate)
Attaching package: ‘lubridate’
The following object is masked from ‘package:cowplot’:
stamp
The following objects are masked from ‘package:Biostrings’:
intersect, setdiff, union
The following objects are masked from ‘package:IRanges’:
%within%, intersect, setdiff, union
The following objects are masked from ‘package:S4Vectors’:
intersect, second, second<-, setdiff, union
The following objects are masked from ‘package:BiocGenerics’:
intersect, setdiff, union
The following objects are masked from ‘package:base’:
date, intersect, setdiff, union
metadata <- read_csv("sample_data.csv")
[36m──[39m [1m[1mColumn specification[1m[22m [36m────────────────────────────────────────────────────────────────────────────────────[39m
cols(
SampleID = [31mcol_character()[39m,
`Year.Trawl#` = [31mcol_character()[39m,
Datecode = [32mcol_double()[39m,
Date = [31mcol_character()[39m,
Month = [32mcol_double()[39m,
Year = [32mcol_double()[39m,
Bayside = [31mcol_character()[39m,
Station = [31mcol_character()[39m
)
Import count table and taxonomy file. I slightly modified otutable.csv in Excel to otutable_mod.csv to remove the quotes around seq names and put NA placehoder as first col name (which was above row names)
# Import Count table. Skip first row of tsv file, which is just some text
count_table <- read_table2("results/otutable_mod.csv")
Missing column names filled in: 'X1' [1]
[36m──[39m [1m[1mColumn specification[1m[22m [36m────────────────────────────────────────────────────────────────────────────────────[39m
cols(
.default = col_double(),
X1 = [31mcol_character()[39m
)
[36mℹ[39m Use [38;5;235m[48;5;253m[38;5;235m[48;5;253m`spec()`[48;5;253m[38;5;235m[49m[39m for the full column specifications.
colnames(count_table)[1] <- "SampleID"
# Import taxonomy of ASVs
taxonomy <- read_csv(file="results/tax_sequences_blast_taxonomy.csv")
Missing column names filled in: 'X1' [1]Duplicated column names deduplicated: 'RefSeq_Tax_ID' => 'RefSeq_Tax_ID_1' [18]
[36m──[39m [1m[1mColumn specification[1m[22m [36m────────────────────────────────────────────────────────────────────────────────────[39m
cols(
X1 = [32mcol_double()[39m,
ASV_ID = [31mcol_character()[39m,
ref_seq_ID = [31mcol_character()[39m,
PID = [32mcol_double()[39m,
alnmt_len = [32mcol_double()[39m,
mismatch = [32mcol_double()[39m,
eval = [32mcol_double()[39m,
bscore = [32mcol_double()[39m,
RefSeq_Tax_ID = [32mcol_double()[39m,
Ref_Seq_title = [31mcol_character()[39m,
superkingdom = [31mcol_character()[39m,
phylum = [31mcol_character()[39m,
class = [31mcol_character()[39m,
order = [31mcol_character()[39m,
family = [31mcol_character()[39m,
genus = [31mcol_character()[39m,
species = [31mcol_character()[39m,
RefSeq_Tax_ID_1 = [32mcol_double()[39m
)
# remove first col of sequential numbers
taxonomy[,1] <- NULL
# filter out sequences with low PID (recommended by Sara)
taxonomy <- filter(taxonomy, PID > 92)
# remove BLAST metadata and just retain taxonomy (necessary for further processing below)
drop.cols <- c(colnames(taxonomy)[2:9],'RefSeq_Tax_ID_1')
taxonomy <- select(taxonomy, -one_of(drop.cols))
# And import the Common names, as curated by Sara. Join to taxonomy
commonnames <- read_excel("Trawls MASTER 2020 _mod_ES.xlsx",7)
commonnames
taxonomy <- left_join(taxonomy, commonnames, by = "ASV_ID")
taxonomy
NA
Filtering removed seqs 110, 332 (Gobiosoma ginsburgi and Belone belone) Note for Sara should we consider setting this at 97% which is more robust and still leaves 334 unique ASVs (rather than 379 with the 92% cutoff in the settings above)
Preview datasets
count_table
taxonomy
metadata
I want to use the phyloseq package for some plotting/ statistics, which first requires making phyloseq objects out of each of input data tables-
count_table_matrix <- as.matrix(count_table[,2:392]) # convert count table to matrix, leaving out character column of sample ID
rownames(count_table_matrix) <- count_table$SampleID # add back in Sample IDs as row names
ASV = otu_table(count_table_matrix, taxa_are_rows = FALSE)
taxonomy_matrix <- as.matrix(taxonomy[,2:9])
rownames(taxonomy_matrix) <- taxonomy$ASV_ID
TAX = tax_table(taxonomy_matrix)
META = sample_data(data.frame(metadata, row.names = metadata$`SampleID`))
First check that the inputs are in compatible formats by checking for ASV names with the phyloseq function, taxa_names
head(taxa_names(TAX))
[1] "Seq_1" "Seq_2" "Seq_3" "Seq_4" "Seq_5" "Seq_6"
head(taxa_names(ASV))
[1] "Seq_1" "Seq_2" "Seq_3" "Seq_4" "Seq_5" "Seq_6"
And check sample names were also detected
# Modify taxa names in ASV, which are formatted with the sample ID, underscor, fastq ID. Don't need this fastq ID anymore and want it to match the sample names from metadata
sample_names(ASV) <- sample_names(ASV) %>%
str_replace_all(pattern = "_S[:digit:]+",replacement = "")
head(sample_names(ASV))
[1] "T1PosCon" "T1S10" "T1S11" "T1S1" "T1S2" "T1S3"
head(sample_names(META))
[1] "T1PosCon" "T1S1" "T1S2" "T1S3" "T1S5" "T1S6"
And make the phyloseq object
ps <- phyloseq(ASV, TAX, META)
rarecurve(otu_table(ps), step=50, cex=0.5)
empty rows removed
# save as .eps
setEPS()
postscript("Figures/rarefaction.eps")
rarecurve(otu_table(ps), step=50, cex=0.5)
empty rows removed
dev.off()
quartz_off_screen
2
Most samples look like they were sampled to completion. Be weary of T3S11, T1S2, and maybe T4S5
Check some features of the phyloseq object
rank_names(ps)
[1] "superkingdom" "phylum" "class" "order" "family" "genus" "species"
[8] "CommonName"
unique(tax_table(ps)[, "superkingdom"])
Taxonomy Table: [2 taxa by 1 taxonomic ranks]:
superkingdom
Seq_1 "Eukaryota"
Seq_377 NA
unique(tax_table(ps)[, "phylum"])
Taxonomy Table: [3 taxa by 1 taxonomic ranks]:
phylum
Seq_1 "Chordata"
Seq_368 "Arthropoda"
Seq_377 NA
unique(tax_table(ps)[, "class"])
Taxonomy Table: [5 taxa by 1 taxonomic ranks]:
class
Seq_1 "Actinopteri"
Seq_63 "Mammalia"
Seq_362 "Chondrichthyes"
Seq_368 "Insecta"
Seq_377 NA
There are some ASVs with NA as superkingdom, phylum, or class annotation- delete these.
ps <- subset_taxa(ps, !is.na(superkingdom) & !is.na(phylum) & !is.na(class))
unique(tax_table(ps)[, "superkingdom"])
Taxonomy Table: [1 taxa by 1 taxonomic ranks]:
superkingdom
Seq_1 "Eukaryota"
unique(tax_table(ps)[, "phylum"])
Taxonomy Table: [2 taxa by 1 taxonomic ranks]:
phylum
Seq_1 "Chordata"
Seq_368 "Arthropoda"
unique(tax_table(ps)[, "class"])
Taxonomy Table: [4 taxa by 1 taxonomic ranks]:
class
Seq_1 "Actinopteri"
Seq_63 "Mammalia"
Seq_362 "Chondrichthyes"
Seq_368 "Insecta"
nrow(tax_table(ps)) # number of ASVs left
[1] 378
378 ASVs still remain…
Also check class Mammalia, to see if contamination or real:
tax_table(subset_taxa(ps, class == 'Mammalia'))
Taxonomy Table: [8 taxa by 8 taxonomic ranks]:
superkingdom phylum class order family genus species CommonName
Seq_63 "Eukaryota" "Chordata" "Mammalia" "Primates" "Hominidae" "Homo" "Homo sapiens" "Human"
Seq_88 "Eukaryota" "Chordata" "Mammalia" "Artiodactyla" "Suidae" "Sus" "Sus scrofa" "Wild boar"
Seq_157 "Eukaryota" "Chordata" "Mammalia" "Primates" "Hominidae" "Homo" "Homo sapiens" "Human"
Seq_343 "Eukaryota" "Chordata" "Mammalia" "Carnivora" "Felidae" "Felis" "Felis catus" "Cat"
Seq_369 "Eukaryota" "Chordata" "Mammalia" "Artiodactyla" "Bovidae" "Bos" "Bos taurus" "Cattle"
Seq_378 "Eukaryota" "Chordata" "Mammalia" "Primates" "Hominidae" "Homo" "Homo sapiens" "Human"
Seq_383 "Eukaryota" "Chordata" "Mammalia" "Primates" "Hominidae" "Homo" "Homo sapiens" "Human"
Seq_389 "Eukaryota" "Chordata" "Mammalia" "Primates" "Hominidae" "Homo" "Homo sapiens" "Human"
These are human, wild boar, cat (…cat lady), and cattle. All are contamination so delete all Mammalia
ps <- subset_taxa(ps, !class == 'Mammalia')
unique(tax_table(ps)[, "class"])
Taxonomy Table: [3 taxa by 1 taxonomic ranks]:
class
Seq_1 "Actinopteri"
Seq_362 "Chondrichthyes"
Seq_368 "Insecta"
Next check the “Insecta” entries
tax_table(subset_taxa(ps, class == 'Insecta'))
Taxonomy Table: [2 taxa by 8 taxonomic ranks]:
superkingdom phylum class order family genus species
Seq_368 "Eukaryota" "Arthropoda" "Insecta" "Hymenoptera" "Formicidae" "Linepithema" "Linepithema humile"
Seq_380 "Eukaryota" "Arthropoda" "Insecta" "Hymenoptera" "Formicidae" "Linepithema" "Linepithema humile"
CommonName
Seq_368 "Ant"
Seq_380 "Ant"
The onlly Insecta is Linepithema humile, which are ants so delete these too..
ps <- subset_taxa(ps, !class == 'Insecta')
unique(tax_table(ps)[, "class"])
Taxonomy Table: [2 taxa by 1 taxonomic ranks]:
class
Seq_1 "Actinopteri"
Seq_362 "Chondrichthyes"
Check overall how the phyla are distributed among samples
# First aglomerate the ASVs at the phylum level using the phyloseq function, tax_glom
superkingdomGlommed = tax_glom(ps, "superkingdom")
# and plot
plot_bar(superkingdomGlommed, x = "Sample")
ggsave(filename = "Figures/seqdepth.eps", plot = plot_bar(superkingdomGlommed, x = "Sample"), units = c("in"), width = 9, height = 6, dpi = 300, )# and save
Total sequences reveals certain samples had very low sequencing effort: T1S7, T1S8, T3S11, and, not as bad, T1S2 and T4S5
The rarefaction analysis also showed T1S2 and T4S5 samples were likely not sequenced to completion. Therefore remove these 5 samples from analysis
ps <- subset_samples(ps, !SampleID == "T1S7" & !SampleID == "T1S8" & !SampleID == "T3S11" & !SampleID == "T1S2" & !SampleID == "T4S5")
ps
phyloseq-class experiment-level object
otu_table() OTU Table: [ 368 taxa and 50 samples ]
sample_data() Sample Data: [ 50 samples by 8 sample variables ]
tax_table() Taxonomy Table: [ 368 taxa by 8 taxonomic ranks ]
50 samples remaining with 368 ASVs
Remove Pos Controls (all hits in positive controls are the same family- I assume this is expected)
ps <- subset_samples(ps, !SampleID == "T1PosCon" & !SampleID == "T2PosCon" & !SampleID == "T3PosCon")
ps
phyloseq-class experiment-level object
otu_table() OTU Table: [ 368 taxa and 47 samples ]
sample_data() Sample Data: [ 47 samples by 8 sample variables ]
tax_table() Taxonomy Table: [ 368 taxa by 8 taxonomic ranks ]
And lastly, correct some taxonomy: According to Sara, Engraulis encrasicolus (European anchovy) should be Anchoa mitchilli (Bay anchovy):
tax_table(ps) <- gsub(tax_table(ps), pattern = "Engraulis encrasicolus", replacement = "Anchoa mitchilli")
47 samples remainwith 368 unique ASVs
For plotting, use relative abundances (# of ASV sequences/sum total sequences in sample), calculated easily using microbiome::transform
ps_ra <- microbiome::transform(ps, transform = "compositional")
Export the relative abundance matrix so Sara can have it:
# Extract abundance matrix from the phyloseq object
RelAbun_matrix = as(otu_table(ps_ra), "matrix")
# Coerce to data.frame
RelAbun_dataframe = as.data.frame(RelAbun_matrix)
# Export
write.csv(RelAbun_dataframe,"results/otutable_relabun.csv", row.names = TRUE)
Then aglomerate the ASVs at the family level using the phyloseq function, tax_glom
familyGlommed_RA = tax_glom(ps_ra, "family")
family_barplot <- plot_bar(familyGlommed_RA, x = "Sample", fill = "family")
family_barplot
NOTES for Sara
Agglomerate by species to see if I get the same 38 unique species Sara sees:
speciesGlommed_RA = tax_glom(ps_ra, "CommonName")
speciesGlommed_RA
phyloseq-class experiment-level object
otu_table() OTU Table: [ 43 taxa and 47 samples ]
sample_data() Sample Data: [ 47 samples by 8 sample variables ]
tax_table() Taxonomy Table: [ 43 taxa by 8 taxonomic ranks ]
tax_table(speciesGlommed_RA)
Taxonomy Table: [43 taxa by 8 taxonomic ranks]:
superkingdom phylum class order family
Seq_1 "Eukaryota" "Chordata" "Actinopteri" "Atheriniformes" "Atherinopsidae"
Seq_2 "Eukaryota" "Chordata" "Actinopteri" "Clupeiformes" "Clupeidae"
Seq_3 "Eukaryota" "Chordata" "Actinopteri" "Clupeiformes" "Engraulidae"
Seq_4 "Eukaryota" "Chordata" "Actinopteri" "Scombriformes" "Pomatomidae"
Seq_5 "Eukaryota" "Chordata" "Actinopteri" "Lutjaniformes" "Lutjanidae"
Seq_6 "Eukaryota" "Chordata" "Actinopteri" "Pleuronectiformes" "Paralichthyidae"
Seq_7 "Eukaryota" "Chordata" "Actinopteri" "Clupeiformes" "Clupeidae"
Seq_9 "Eukaryota" "Chordata" "Actinopteri" "Gobiiformes" "Gobiidae"
Seq_10 "Eukaryota" "Chordata" "Actinopteri" "Pleuronectiformes" "Scophthalmidae"
Seq_11 "Eukaryota" "Chordata" "Actinopteri" "Perciformes" "Serranidae"
Seq_12 "Eukaryota" "Chordata" "Actinopteri" "Spariformes" "Sparidae"
Seq_15 "Eukaryota" "Chordata" "Actinopteri" NA "Sciaenidae"
Seq_16 "Eukaryota" "Chordata" "Actinopteri" NA "Sciaenidae"
Seq_17 "Eukaryota" "Chordata" "Actinopteri" "Labriformes" "Labridae"
Seq_19 "Eukaryota" "Chordata" "Actinopteri" "Perciformes" "Cottidae"
Seq_20 "Eukaryota" "Chordata" "Actinopteri" "Pleuronectiformes" "Pleuronectidae"
Seq_21 "Eukaryota" "Chordata" "Actinopteri" NA "Moronidae"
Seq_22 "Eukaryota" "Chordata" "Actinopteri" "Syngnathiformes" "Syngnathidae"
Seq_30 "Eukaryota" "Chordata" "Actinopteri" "Pleuronectiformes" "Paralichthyidae"
Seq_33 "Eukaryota" "Chordata" "Actinopteri" NA "Sciaenidae"
Seq_34 "Eukaryota" "Chordata" "Actinopteri" "Labriformes" "Labridae"
Seq_36 "Eukaryota" "Chordata" "Actinopteri" "Anguilliformes" "Anguillidae"
Seq_38 "Eukaryota" "Chordata" "Actinopteri" "Scombriformes" "Scombridae"
Seq_40 "Eukaryota" "Chordata" "Actinopteri" "Perciformes" "Gasterosteidae"
Seq_44 "Eukaryota" "Chordata" "Actinopteri" "Cyprinodontiformes" "Fundulidae"
Seq_50 "Eukaryota" "Chordata" "Actinopteri" "Atheriniformes" "Atherinopsidae"
Seq_52 "Eukaryota" "Chordata" "Actinopteri" "Gadiformes" "Phycidae"
Seq_54 "Eukaryota" "Chordata" "Actinopteri" "Scombriformes" "Scombridae"
Seq_57 "Eukaryota" "Chordata" "Actinopteri" "Perciformes" "Triglidae"
Seq_67 "Eukaryota" "Chordata" "Actinopteri" "Scombriformes" "Scombridae"
Seq_82 "Eukaryota" "Chordata" "Actinopteri" NA "Sciaenidae"
Seq_84 "Eukaryota" "Chordata" "Actinopteri" "Gadiformes" "Gadidae"
Seq_102 "Eukaryota" "Chordata" "Actinopteri" "Clupeiformes" "Engraulidae"
Seq_103 "Eukaryota" "Chordata" "Actinopteri" "Perciformes" "Cottidae"
Seq_115 "Eukaryota" "Chordata" "Actinopteri" "Cyprinodontiformes" "Fundulidae"
Seq_119 "Eukaryota" "Chordata" "Actinopteri" "Gadiformes" "Phycidae"
Seq_139 "Eukaryota" "Chordata" "Actinopteri" "Batrachoidiformes" "Batrachoididae"
Seq_141 "Eukaryota" "Chordata" "Actinopteri" "Scombriformes" "Scombridae"
Seq_181 "Eukaryota" "Chordata" "Actinopteri" "Tetraodontiformes" "Tetraodontidae"
Seq_231 "Eukaryota" "Chordata" "Actinopteri" "Gadiformes" "Merlucciidae"
Seq_359 "Eukaryota" "Chordata" "Actinopteri" "Perciformes" "Triglidae"
Seq_362 "Eukaryota" "Chordata" "Chondrichthyes" "Myliobatiformes" "Myliobatidae"
Seq_372 "Eukaryota" "Chordata" "Chondrichthyes" "Carcharhiniformes" "Triakidae"
genus species CommonName
Seq_1 "Menidia" "Menidia menidia" "Atlantic silverside"
Seq_2 "Brevoortia" "Brevoortia tyrannus" "Atlantic menhaden"
Seq_3 "Engraulis" "Anchoa mitchilli" "Bay anchovy"
Seq_4 "Pomatomus" "Pomatomus saltatrix" "Bluefish"
Seq_5 "Lutjanus" "Lutjanus griseus" "Grey snapper"
Seq_6 "Paralichthys" "Paralichthys dentatus" "Summer flounder"
Seq_7 "Alosa" "Alosa mediocris" "Hickory shad"
Seq_9 "Gobiosoma" "Gobiosoma ginsburgi" "Seaboard goby"
Seq_10 "Scophthalmus" "Scophthalmus aquosus" "Windowpane flounder"
Seq_11 "Centropristis" "Centropristis striata" "Black seabass"
Seq_12 "Stenotomus" "Stenotomus chrysops" "Scup"
Seq_15 "Leiostomus" "Leiostomus xanthurus" "Spot"
Seq_16 "Menticirrhus" "Menticirrhus saxatilis" "Northern kingfish"
Seq_17 "Tautoga" "Tautoga onitis" "Tautog"
Seq_19 "Myoxocephalus" "Myoxocephalus aenaeus" "Grubby sculpin"
Seq_20 "Pseudopleuronectes" "Pseudopleuronectes americanus" "Winter flounder"
Seq_21 "Morone" "Morone saxatilis" "Striped bass"
Seq_22 "Syngnathus" "Syngnathus fuscus" "Northern pipefish"
Seq_30 "Etropus" "Etropus microstomus" "Smallmouth flounder"
Seq_33 "Cynoscion" "Cynoscion regalis" "Weakfish"
Seq_34 "Tautogolabrus" "Tautogolabrus adspersus" "Cunner"
Seq_36 "Anguilla" "Anguilla rostrata" "American eel"
Seq_38 "Thunnus" "Thunnus obesus" "Bigeye tuna"
Seq_40 "Apeltes" "Apeltes quadracus" "Stickleback"
Seq_44 "Fundulus" "Fundulus majalis" "Striped killifish"
Seq_50 "Membras" "Membras martinica" "Rough silverside"
Seq_52 "Urophycis" "Urophycis floridana" "Spotted hake"
Seq_54 "Scomber" "Scomber japonicus" "Chub mackerel"
Seq_57 "Prionotus" "Prionotus carolinus" "Northern searobin"
Seq_67 "Thunnus" "Thunnus thynnus" "Atlantic bluefin tuna"
Seq_82 "Bairdiella" "Bairdiella chrysoura" "American silver perch"
Seq_84 "Microgadus" "Microgadus tomcod" "Atlantic tomcod"
Seq_102 "Engraulis" "Engraulis mordax" "Bay anchovy"
Seq_103 "Myoxocephalus" "Myoxocephalus quadricornis" "Fourhorn sculpin"
Seq_115 "Fundulus" "Fundulus heteroclitus" "Mummichog"
Seq_119 "Urophycis" "Urophycis floridana" "Red hake"
Seq_139 "Opsanus" "Opsanus tau" "Oyster toadfish"
Seq_141 "Katsuwonus" "Katsuwonus pelamis" "Skipjack tuna"
Seq_181 "Sphoeroides" "Sphoeroides maculatus" "Northern puffer"
Seq_231 "Merluccius" "Merluccius bilinearis" "Silver hake"
Seq_359 "Prionotus" "Prionotus evolans" "Striped searobin"
Seq_362 "Rhinoptera" "Rhinoptera bonasus" "Cownose ray"
Seq_372 "Mustelus" "Mustelus canis" "Dusky smooth-hound shark"
NOTES for Sara
Based on my previous scripts with Cariaco Eukaryotic data
# convert ps object to dataframe using phyloseq's psmelt
species_df <- psmelt(speciesGlommed_RA)
# replace zeroes in the table with NA
species_df[species_df == 0] <- NA
# and remove rows with NAs in abundance (this is so they don't appear as small dots in plot)
species_df <- filter(species_df, !is.na(Abundance))
Plot by species, scientific name
speciesbubbleplot_eDNA_sciname <- ggplot(species_df, aes(x = Station, y = fct_rev(species), color = Station)) + # the fancy stuff around y (species) helps to present it in reverse order in the plot (from top to btm alphabetically)
geom_point(aes(size = Abundance, fill = Station), color = "black", pch = 21)+
scale_size(range = c(1,15)) +
scale_size_area(breaks = c(0,.25,.5,.75,1), max_size = 6)+
xlab("")+
ylab("")+
labs(size="Relative Abundance")+
theme_bw() +
scale_fill_brewer(palette="Paired") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
facet_grid(Datecode~Bayside, scales = "free", space = "free", drop= TRUE)
Scale for 'size' is already present. Adding another scale for 'size', which will replace the existing
scale.
speciesbubbleplot_eDNA_sciname
Plot by species common name
speciesbubbleplot_eDNA_comname <- ggplot(species_df, aes(x = Station, y = fct_rev(CommonName), color = Station)) + # the fancy stuff around y (CommonName) helps to present it in reverse order in the plot (from top to btm alphabetically)
geom_point(aes(size = Abundance, fill = Station), color = "black", pch = 21)+
scale_size(range = c(1,15)) +
scale_size_area(breaks = c(0,.25,.5,.75,1), max_size = 6)+
xlab("")+
ylab("")+
labs(size="Relative Abundance")+
theme_bw() +
scale_fill_brewer(palette="Paired") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
facet_grid(Datecode~Bayside, scales = "free", space = "free", drop= TRUE)
Scale for 'size' is already present. Adding another scale for 'size', which will replace the existing
scale.
speciesbubbleplot_eDNA_comname
Exportfigures
ggsave(filename = "Figures/speciesbubbleplot_eDNA_sciname.eps", plot = speciesbubbleplot_eDNA_sciname, units = c("in"), width = 7, height = 12, dpi = 300)
ggsave(filename = "Figures/speciesbubbleplot_eDNA_comname.eps", plot = speciesbubbleplot_eDNA_comname, units = c("in"), width = 7, height = 12, dpi = 300)
NOTE on above. The common name plot has two entries in the Bay anchovy row because, as mentioned above, there are two different species name that are labelled as Bay Anchovy. Is it OK to group these as same species (Anchoa mitchilli)
# import 4th sheet from Excel file which contains morphometric data for each individual collected for every date
trawl_master <- read_excel("Trawls MASTER 2020 _mod_ES.xlsx",4)
trawl_master
# and import 6th sheet which is station info
stations <- read_excel("Trawls MASTER 2020 _mod_ES.xlsx",6)
stations
NA
Make an equivalent to an OTU table, grouping by date and location and representing counts for every unique species
trawl_counts <- trawl_master %>%
group_by(DATECODE, STATION_NO, CommonName) %>%
tally(name = "count")
trawl_counts
and link station names instead of numbers to count table
trawl_counts <- left_join(trawl_counts, stations, by = "STATION_NO")
trawl_counts
Remove 09/16/20 since there is no equivalent eDNA from that date
trawl_counts <- trawl_counts %>%
filter(DATECODE != "20200916")
trawl_counts
speciesbubbleplot_trawl_comname <- ggplot(trawl_counts, aes(x = STATION_NA, y = fct_rev(CommonName), color = STATION_NA)) +
geom_point(aes(size = log10(count), fill = STATION_NA), color = "black", pch = 21)+
scale_size(range = c(1,15)) +
scale_size_area(breaks = c(.01,.1, .3, .5, 1, 3), max_size = 6)+
xlab("")+
ylab("")+
labs(size="Log(counts)", fill = "Station")+
theme_bw() +
scale_fill_brewer(palette="Paired") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
facet_grid(DATECODE~BAYSIDE, scales = "free", space = "free", drop= TRUE)
Scale for 'size' is already present. Adding another scale for 'size', which will replace the existing
scale.
speciesbubbleplot_trawl_comname
Export figure
ggsave(filename = "Figures/speciesbubbleplot_trawl_comname.eps", plot = speciesbubbleplot_trawl_comname, units = c("in"), width = 6.75, height = 13, dpi = 300)
Count unique species across all stations, grouped by date, for each method (trawl, eDNA)
trawl_uniques <- trawl_counts %>%
group_by(DATECODE, CommonName) %>%
summarise(Trawl_Count = sum(count, na.rm=TRUE))
`summarise()` regrouping output by 'DATECODE' (override with `.groups` argument)
trawl_uniques
eDNA_uniques <- species_df%>%
group_by(Datecode, CommonName) %>%
summarise(eDNA_RelAbun = sum(Abundance, na.rm=TRUE))
`summarise()` regrouping output by 'Datecode' (override with `.groups` argument)
eDNA_uniques
# Combine into one dataframe
trawl_eDNA_abun_table <- full_join(trawl_uniques, eDNA_uniques, by=c("CommonName" = "CommonName", "DATECODE" = "Datecode"))
trawl_eDNA_abun_table
Count total number of species from each method for each date
eDNA_richness <- tally(eDNA_uniques, name = "eDNA")
trawl_richness <- tally(trawl_uniques, name = "trawl")
speciesrichness <- full_join(eDNA_richness, trawl_richness, c("Datecode" = "DATECODE"))
speciesrichness <- pivot_longer(speciesrichness, !Datecode, names_to = "Method", values_to = "Richness")
speciesrichness$Datecode <- ymd(speciesrichness$Datecode) # convert to date format (better for plotting)
speciesrichness
Plot side-by-side
species_richness_plot <- ggplot(speciesrichness, aes(x =Datecode, y = Richness)) +
geom_line(aes(color = Method), size = 3) +
theme_bw() +
xlab("") +
ylab("Species Richness")
species_richness_plot
# export plot
ggsave(filename = "Figures/species_richness_plot.eps", plot = species_richness_plot, units = c("in"), width = 4, height = 3, dpi = 300)
Sum total number of species across all dates/ stations for entire study
species_sums_abun_table <- trawl_eDNA_abun_table %>%
group_by(CommonName) %>%
summarise(Trawl = sum(Trawl_Count, na.rm=TRUE), eDNA = (sum(eDNA_RelAbun, na.rm=TRUE))) %>%
pivot_longer(!CommonName, names_to = "Method", values_to = "Abundance")
`summarise()` ungrouping output (override with `.groups` argument)
# turn zeroes to NA so they don't plot
species_sums_abun_table <- na_if(species_sums_abun_table,0)
species_sums_abun_table
For each species, plot side-by-side comparison of abundance (summed over whole study) using each method
# First create a custom color scale to make this pretty
myColors <- colorRampPalette(brewer.pal(11,"Spectral"))(55)
names(myColors) <- levels(unique(species_sums_abun_table$CommonName))
colScale <- scale_colour_manual(name = "CommonName",values = myColors)
species_abun_sum_plot <- ggplot(species_sums_abun_table, aes(x = Abundance, y = reorder(CommonName, Abundance, function(x){sum(x,na.rm = TRUE)}), color = CommonName)) +
geom_point(size = 5) +
facet_wrap(~fct_rev(Method), scales = "free") +
theme_bw() +
xlab("Abundance") +
ylab("") +
colScale +
theme(legend.position = "none")
species_abun_sum_plot
Export plot
ggsave(filename = "Figures/species_abun_sum_plot.eps", plot = species_abun_sum_plot, units = c("in"), width = 7, height = 8, dpi = 300)
PCA is essentially a type of PCoA using the Euclidean distance matrix as input. When combined with a log-ratio transformation of the count table, this is deemed appropriate for compositional datasets. It is also recommended as a first step in exploratory analyses of sequencinging datasets.
First do a CLR, centered log ratio transformation of the absolute abundance data (after filtering), as suggested by Gloor et al. 2017
# Estimate covariance matrix for CLR-transformed ASV table
clr_asv_table_ps <- data.frame(compositions::clr(otu_table(ps)))
Generate the PCA and visualize axes
# Generate a Principle Component Analysis (PCA) and evaluated based on the eigen decomposition from sample covariance matrix.
lograt_pca <- prcomp(clr_asv_table_ps)
# NOTE- this is equivalent to first making a Euclidean distance matrix using the CLR data table and then running a PCoA. A Euclidean distance matrix of a log-transformed data table = an Aitchison distance matrix. So this is equivalent to the compositional methods listed in Gloor et al.
# Visual representation with a screeplot
lograt_variances <- as.data.frame(lograt_pca$sdev^2/sum(lograt_pca$sdev^2)) %>% #Extract axes
# Format to plot
select(PercVar = 'lograt_pca$sdev^2/sum(lograt_pca$sdev^2)') %>%
rownames_to_column(var = "PCaxis") %>%
data.frame
head(lograt_variances)
# Plot screeplot
ggplot(lograt_variances, aes(x = as.numeric(PCaxis), y = PercVar)) +
geom_bar(stat = "identity", fill = "grey", color = "black") +
theme_minimal() +
theme(axis.title = element_text(color = "black", face = "bold", size = 10),
axis.text.y = element_text(color = "black", face = "bold"),
axis.text.x = element_blank()) +
labs(x = "PC axis", y = "% Variance", title = "Log-Ratio PCA Screeplot, CLR Tranformation")
Plot in 3D using first 3 axes since the 2nd and 3rd are similar proportions of variance. Total variance explained by first three: 15.7 + 10.5 + 10.0 = 36.2%)
Visualize the PCA-
# Extract variances from the clr pca
There were 14 warnings (use warnings() to see them)
pca_lograt_frame <- data.frame(lograt_pca$x) %>%
rownames_to_column(var = "SampleID")
# Merge metadata into the pcoa data table
pca_lograt_frame <- left_join(pca_lograt_frame, metadata, by = "SampleID")
head(pca_lograt_frame)
# Select eigenvalues from dataframe, round to 4 places and multiply by 100 for plotting. These will be the axes for the 3-D plot
eigenvalues<-round(lograt_variances[,2], digits = 4)*100
# Plotly - 3-D
pca_lograt <- plot_ly(pca_lograt_frame, type='scatter3d', mode='markers',
x=~PC1,y=~PC2,z=~PC3,colors=~brewer.pal(11,'Paired'),
color=~Station, symbols = c('circle','diamond'), symbol=~Bayside)%>%
layout(font=list(size=12),
title='CLR-Euclidean PCA',
scene=list(xaxis=list(title=paste0('Co 2 ',eigenvalues[2],'%'),
showticklabels=FALSE,zerolinecolor='black'),
yaxis=list(title=paste0('Co 3 ',eigenvalues[3],'%'),
showticklabels=FALSE,zerolinecolor='black'),
zaxis=list(title=paste0('Co 1 ',eigenvalues[1],'%'),
showticklabels=FALSE,zerolinecolor='black')))
pca_lograt
# save in "Embedded_figures" dirctory so that it can be hosted at Github and embedded in this notebook
withr::with_dir('Embedded_figures', htmlwidgets::saveWidget(as_widget(pca_lograt), file="pca_lograt.html", selfcontained = F))
# plot so that it renders in github notebook
htmltools::tags$iframe(
src=file.path("Embedded_figures/pca_lograt.html"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0"
)
NA
The CLR-Euclidean PCA reveals there is some separation according to East vs West. The PCA only explains ~36% of the variance so keep going with different ordinations to see if we can get a better representation
Test embedding image in-line instead of as code chunk:
The more traditional approach to ordinations is to do a PCoA on a distance matrix such as Bray-Curtis, Jaccard, or Unifrac. While these are not considered compositional approaches, when combined with pre-treatment (transformations) they become more appropriate. One such transformation that I will use here is the Hellinger transformation.
The different distance matrices also tell you a few different things about the dataset so I will run through this to try to see if I can tease those out.
Before calculating any distance matrix, do a transformation of the filtered count table. Hellinger transformation is the square root of the relative abundance, so calculate it based on the ps_ra object:
ps_hellinger <- transform_sample_counts(ps_ra, function(x){sqrt(x)})
First, Jaccard, which builds the distance matrix based on presence/absence between samples. It does not take into account relative abundance of the taxa. Therefore this functions well for determining differences driven by rare taxa, which are weighed the same as abundant taxa.
jac_dmat<-vegdist(otu_table(ps_hellinger),method="jaccard") # Jaccard dist metric
pcoa_jac<-ape::pcoa(jac_dmat) # perform PCoA
# Extract variances from pcoa, from jaccard calculated dist. metric
jac_variances <- data.frame(pcoa_jac$values$Relative_eig) %>%
select(PercVar = 'pcoa_jac.values.Relative_eig') %>%
rownames_to_column(var = "PCaxis") %>%
data.frame
head(jac_variances)
# Make a screeplot
ggplot(jac_variances, aes(x = as.numeric(PCaxis), y = PercVar)) +
geom_bar(stat = "identity", fill = "grey", color = "black") +
theme_minimal() +
theme(axis.title = element_text(color = "black", face = "bold", size = 10),
axis.text.y = element_text(color = "black", face = "bold"),
axis.text.x = element_blank()) +
labs(x = "PC axis", y = "% Variance", title = "Jaccard PCoA Screeplot")
The first two axes (19.0 + 9.6 = 28.6) are OK. But I am going to experiment and plot the first 3 axes since the 2nd and 3rd explain a similar amount of variance, 19.6 and 8.4% respectively
Plot in 3D with Plotly
# Extract variances from the jaccard pcoa
pcoa_jac_df <- data.frame(pcoa_jac$vectors) %>%
rownames_to_column(var = "SampleID")
# Merge metadata into the pcoa data table
pcoa_jac_df <- left_join(pcoa_jac_df, metadata, by = "SampleID")
head(pcoa_jac_df)
# Select eigenvalues from dataframe, round to 4 places and multiply by 100 for plotting. These will be the axes for the 3-D plot
eigenvalues<-round(jac_variances[,2], digits = 4)*100
# Plotly - 3-D
pcoa_jaccard <- plot_ly(pcoa_jac_df, type='scatter3d', mode='markers',
x=~Axis.2,y=~Axis.3,z=~Axis.1,colors=~brewer.pal(11,'Paired'),
color=~Station, symbols = c('circle','diamond'), symbol=~Bayside)%>%
layout(font=list(size=12),
title='PCoA Jaccard Distance',
scene=list(xaxis=list(title=paste0('Co 2 ',eigenvalues[2],'%'),
showticklabels=FALSE,zerolinecolor='black'),
yaxis=list(title=paste0('Co 3 ',eigenvalues[3],'%'),
showticklabels=FALSE,zerolinecolor='black'),
zaxis=list(title=paste0('Co 1 ',eigenvalues[1],'%'),
showticklabels=FALSE,zerolinecolor='black')))
pcoa_jaccard
`arrange_()` is deprecated as of dplyr 0.7.0.
Please use `arrange()` instead.
See vignette('programming') for more help
[90mThis warning is displayed once every 8 hours.[39m
[90mCall `lifecycle::last_warnings()` to see where this warning was generated.[39m
withr::with_dir('Figures', htmlwidgets::saveWidget(as_widget(pcoa_jaccard), file="pcoa_jaccard.html"))
The Jaccard-PCoA shows separation along axis 2 in East vs West differences.
Next, try a Bray-Curtis distance matrix with PCoA, which builds the distance matrix based on presence/absence between samples and relative abundance differences. This ordination will represent well the differences in samples that are driven by taxa with high relative abundances.
bray_dmat<-vegdist(otu_table(ps_hellinger),method="bray") # Bray-Curtis dist metric
pcoa_bray<-ape::pcoa(bray_dmat) # perform PCoA
# Extract variances from pcoa, from jaccard calculated dist. metric
bray_variances <- data.frame(pcoa_bray$values$Relative_eig) %>%
select(PercVar = 'pcoa_bray.values.Relative_eig') %>%
rownames_to_column(var = "PCaxis") %>%
data.frame
head(bray_variances)
# Make a screeplot
ggplot(bray_variances, aes(x = as.numeric(PCaxis), y = PercVar)) +
geom_bar(stat = "identity", fill = "grey", color = "black") +
theme_minimal() +
theme(axis.title = element_text(color = "black", face = "bold", size = 10),
axis.text.y = element_text(color = "black", face = "bold"),
axis.text.x = element_blank()) +
labs(x = "PC axis", y = "% Variance", title = "Bray-Curtis PCoA Screeplot")
The first two axes (27.7 + 14.3 = 42%) are pretty good again but I am still going to experiment in the plot with the 3rd axis since it is similar to the second (12.2% variance)
Plot in 3D with Plotly
# Extract variances from the jaccard pcoa
pcoa_bray_df <- data.frame(pcoa_bray$vectors) %>%
rownames_to_column(var = "SampleID")
# Merge metadata into the pcoa data table
pcoa_bray_df <- left_join(pcoa_bray_df, metadata, by = "SampleID")
head(pcoa_bray_df)
# Select eigenvalues from dataframe, round to 4 places and multiply by 100 for plotting. These will be the axes for the 3-D plot
eigenvalues<-round(bray_variances[,2], digits = 4)*100
# Plotly - 3-D
pcoa_bray <- plot_ly(pcoa_bray_df, type='scatter3d', mode='markers',
x=~Axis.2,y=~Axis.3,z=~Axis.1,colors=~brewer.pal(11,'Paired'),
color=~Station, symbols = c('circle','diamond'), symbol=~Bayside)%>%
layout(font=list(size=12),
title='PCoA Bray-Curtis Distance',
scene=list(xaxis=list(title=paste0('Co 2 ',eigenvalues[2],'%'),
showticklabels=FALSE,zerolinecolor='black'),
yaxis=list(title=paste0('Co 3 ',eigenvalues[3],'%'),
showticklabels=FALSE,zerolinecolor='black'),
zaxis=list(title=paste0('Co 1 ',eigenvalues[1],'%'),
showticklabels=FALSE,zerolinecolor='black')))
pcoa_bray
withr::with_dir('Figures', htmlwidgets::saveWidget(as_widget(pcoa_bray), file="pcoa_bray.html"))
These results are similar to Jaccard: the second axis seems driven by differences in East vs West. But there are clearly other things going on here with axies 1 and 3. I think this is a good representation of the data: together the 3 axes explain 54.13% of the variance.
Lastly, try a non-metric dimensional scaling ordination. PCA/PCoA are metric and attempt to rotate axes to fit the distance matrix distribution. An NMDS represents the data in 2-axes, by constraining the distribution of the points. Similar to above, this can be combined with different pre-treatment of the data.
First try the compositional approach, an NMDS on CLR-tranformed data using the Euclidean distances (aka Aitchison distance)
euc_dmat<-dist(clr_asv_table_ps, method = "euclidean") # Build the Aitchison distance matrix
euc_nmds <- metaMDS(euc_dmat, k=2, autotransform=FALSE) # Run the ordination
Run 0 stress 0.2105436
Run 1 stress 0.2110814
Run 2 stress 0.2130335
Run 3 stress 0.2337817
Run 4 stress 0.2311286
Run 5 stress 0.2110906
Run 6 stress 0.2224957
Run 7 stress 0.2158067
Run 8 stress 0.2110925
Run 9 stress 0.2233104
Run 10 stress 0.2316347
Run 11 stress 0.2107091
... Procrustes: rmse 0.005974315 max resid 0.03264245
Run 12 stress 0.2314
Run 13 stress 0.2287184
Run 14 stress 0.2127398
Run 15 stress 0.2272337
Run 16 stress 0.2160322
Run 17 stress 0.2320417
Run 18 stress 0.2145287
Run 19 stress 0.2156851
Run 20 stress 0.2110731
*** No convergence -- monoMDS stopping criteria:
1: no. of iterations >= maxit
19: stress ratio > sratmax
euc_nmds$stress #Check the stress. Less than 0.1 is good. Less than 0.05 is better. This will be different each time, since it is iteratively finding a unique solution each time (although the should look similar)
[1] 0.2105436
# Extract points from nmds and merge into data frame with metadata
euc_nmds_df <- data.frame(euc_nmds$points) %>%
rownames_to_column(var = "SampleID")
# Merge metadata into the pcoa data table
euc_nmds_df <- left_join(euc_nmds_df, metadata, by = "SampleID")
head(euc_nmds_df)
## Plotting euclidean distance NMDS
nmds_aitch <- ggplot(euc_nmds_df,aes(x = MDS1, y = MDS2, color = Station, shape = Bayside)) +
geom_point(size = 4) +
scale_color_brewer(palette="Paired") +
theme_bw() +
labs(x = "NMDS 1", y = "NMDS 2", title = paste0('Aitchison Distance NMDS, Stress = ', round(euc_nmds$stress,2))) +
coord_fixed(ratio = 1)
nmds_aitch
ggsave("figures/nmds_aitch.eps",nmds_aitch, width = 7, height = 5, units = c("in"))
The above has a relatively high stress (>0.2) so should be interpreted with caution. But it does show some separation East vs West along NMDS 1.
Next try a Jaccard NMDS, which will represent differences in presence/absence among samples, emphasizing both abundant and rare taxa the same
jac_nmds <- metaMDS(jac_dmat, k=2, autotransform=FALSE) # Run the ordination. Distance matrix was already calculated above
Run 0 stress 0.1627003
Run 1 stress 0.1633416
Run 2 stress 0.1662452
Run 3 stress 0.1625701
... New best solution
... Procrustes: rmse 0.004322226 max resid 0.02355899
Run 4 stress 0.1496845
... New best solution
... Procrustes: rmse 0.09101825 max resid 0.3182899
Run 5 stress 0.1511997
Run 6 stress 0.1573356
Run 7 stress 0.1573408
Run 8 stress 0.1744426
Run 9 stress 0.1574038
Run 10 stress 0.1574293
Run 11 stress 0.1625698
Run 12 stress 0.1713826
Run 13 stress 0.1498312
... Procrustes: rmse 0.01157208 max resid 0.06868134
Run 14 stress 0.150871
Run 15 stress 0.1707482
Run 16 stress 0.1705187
Run 17 stress 0.1508708
Run 18 stress 0.1721511
Run 19 stress 0.1508709
Run 20 stress 0.1496849
... Procrustes: rmse 0.000227522 max resid 0.001291969
... Similar to previous best
*** Solution reached
jac_nmds$stress #Check the stress. Less than 0.1 is good. Less than 0.5 is better. This will be different each time, since it is iteratively finding a unique solution each time (although the should look similar)
[1] 0.1496845
# Extract points from nmds and merge into data frame with metadata
jac_nmds_df <- data.frame(jac_nmds$points) %>%
rownames_to_column(var = "SampleID")
# Merge metadata into the pcoa data table
jac_nmds_df <- left_join(jac_nmds_df, metadata, by = "SampleID")
head(jac_nmds_df)
## Plotting euclidean distance NMDS
nmds_jaccard <- ggplot(jac_nmds_df,aes(x = MDS1, y = MDS2, color = Station, shape = Bayside)) +
geom_point(size = 4) +
scale_color_brewer(palette="Paired") +
theme_bw() +
labs(x = "NMDS 1", y = "NMDS 2", title = paste0('Jaccard Distance NMDS, Stress = ', round(jac_nmds$stress,2))) +
coord_fixed(ratio = 1)
nmds_jaccard
ggsave("figures/nmds_jaccard.eps",nmds_jaccard, width = 7, height = 5, units = c("in"))
This is still a relatively high stress (>0.1) so should be interpreted with caution. Similar to Aitchison-distance nMDS but there is a little more separation of East vs West on NMDS 2 axis.
Next try a Bray-Curis NMDS, which will represent differences in presence/absence among samples and relative abundance, thus emphasizing impacts of highly abundant taxa.
bray_nmds <- metaMDS(bray_dmat, k=2, autotransform=FALSE) # Run the ordination. Distance matrix was already calculated above
Run 0 stress 0.1628608
Run 1 stress 0.1700644
Run 2 stress 0.1706993
Run 3 stress 0.1574041
... New best solution
... Procrustes: rmse 0.08465261 max resid 0.3069557
Run 4 stress 0.1713996
Run 5 stress 0.1805603
Run 6 stress 0.1572102
... New best solution
... Procrustes: rmse 0.02416992 max resid 0.1186635
Run 7 stress 0.1498307
... New best solution
... Procrustes: rmse 0.05787324 max resid 0.3389741
Run 8 stress 0.1764597
Run 9 stress 0.1512066
Run 10 stress 0.1495162
... New best solution
... Procrustes: rmse 0.05231193 max resid 0.3245111
Run 11 stress 0.1652534
Run 12 stress 0.1627925
Run 13 stress 0.1574076
Run 14 stress 0.1495156
... New best solution
... Procrustes: rmse 0.0005084261 max resid 0.003155013
... Similar to previous best
Run 15 stress 0.1498312
... Procrustes: rmse 0.05220017 max resid 0.3259805
Run 16 stress 0.1568971
Run 17 stress 0.1496845
... Procrustes: rmse 0.05079648 max resid 0.3274102
Run 18 stress 0.2011478
Run 19 stress 0.166024
Run 20 stress 0.1634199
*** Solution reached
bray_nmds$stress #Check the stress. Less than 0.1 is good. Less than 0.5 is better. This will be different each time, since it is iteratively finding a unique solution each time (although the should look similar)
[1] 0.1495156
# Extract points from nmds and merge into data frame with metadata
bray_nmds_df <- data.frame(bray_nmds$points) %>%
rownames_to_column(var = "SampleID")
# Merge metadata into the pcoa data table
bray_nmds_df <- left_join(bray_nmds_df, metadata, by = "SampleID")
head(bray_nmds_df)
## Plotting euclidean distance NMDS
nmds_bray <- ggplot(bray_nmds_df,aes(x = MDS1, y = MDS2, color = Station, shape = Bayside)) +
geom_point(size = 4) +
scale_color_brewer(palette="Paired") +
theme_bw() +
labs(x = "NMDS 1", y = "NMDS 2", title = paste0('Bray-Curtis Distance NMDS, Stress = ', round(bray_nmds$stress,2))) +
coord_fixed(ratio = 1)
nmds_bray
ggsave("figures/nmds_bray.eps",nmds_bray, width = 7, height = 5, units = c("in"))
Very similar to Jaccard results. High-ish stress (0.15)
–> CONTINUE HERE
XXX